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ABSTRACT 

We present the two- and three-point real space correlation functions of the five-year WMAP sky maps 
and compare the observed functions to simulated ACDM concordance model ensembles. In agreement 
with previously published results, we find that the temperature correlation functions are consistent 
with expectations. However, the pure polarization correlation functions are acceptable only for the 
33GHz band map; the 41, 61, and 94 GHz band correlation functions all exhibit significant large-scale 
excess structures. Further, these excess structures very closely match the correlation functions of 
the two (synchrotron and dust) foreground templates used to correct the WMAP data for galactic 
contamination, with a cross-correlation statistically significant at the 2a — 3er confidence level. The 
correlation is slightly stronger with respect to the thermal dust template than with the synchrotron 
template. 

Subject headings: cosmic background radiation — cosmology: observations — methods: numerical 



1. INTRODUCTION 

Observations of the cosmic microwave background 
(CMB) have been among the most important ingredi- 
ents in the revolution of cosmology that has taken place 
in the last two decades, when cosmology changed from 
a data starved to a data rich science. The so far most 
influential observations have been made with the Wilkin- 
son Microwave Anisotropy Probe ( WMAP; Bennett et al. 
2003a; Hinshaw et al. 2007, 2009) satellite experiment, 
which has measured the microwave sky at five different 
frequencies (23, 33, 41, 61, and 94 GHz) in both tem- 
perature and polarization. Its main results are five sky 
maps with resolutions between 13' and 55', each compris- 
ing more than three million pixels in each of the three 
Stokes parameters I, Q, and U. 

The detector technology leading to instruments like 
WMAP has thus been very important in observing the 
CMB. Almost equally important has been the progress 
in computer technology and algorithms. With the 
enormously increased number of pixels from new high- 
resolution instruments, scientists have had to develop 
clever algorithms for every step of the required data anal- 
ysis pipeline: map making (e.g., Ashdown et al. 2007; 
Hinshaw et al. 2003b, and references therein), compo- 
nent separation (e.g., Bennett et al. 2003b; Leach et al. 
2008; Eriksen et al. 2008, and references therein), power 
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spectrum and cosmological parameter estimation (e.g., 
Gorski 1994; Lewis & Bridle 2002; Hivon et al. 2002; 
Hinshaw et al. 2003a; Verde et al. 2003; Eriksen et al. 
2004d, and references therein), and analysis of higher- 
order statistics (e.g., Hinshaw et al. 1994; Kogut et al. 
1995; Komatsu et al. 2003; Eriksen et al. 2004c, 2005 
and references therein). 

In this paper, we revisit a well-known example of the 
latter category, namely, real-space TV-point correlation 
functions. These functions arise naturally in studies of 
non-Gaussianity, since it can be shown that any odd- 
ordered TV-point function has an identically vanishing 
expectation value, while all even-ordered TV-point func- 
tions have an expectation value given by products of the 
corresponding two-point functions. Violation of either 
of these relations would indicate the presence of a non- 
Gaussian component in the field under consideration. 

Earlier TV-point correlation function analyses of 
CMB data inc l ude s tudies of the COB E-BMR data 
dHinshaw et all 119951: iKogut et all 119961: 
2002), the WMAP temperature data 
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|2005j) and one single two- point analysis of th e first-year 
WMAP polarization data (|Kogut et all 120031 ). This pa- 
per is the first to consider the much more mature five- 
year WMAP polarization data, and the first to compute 
the three-point correlation function from any CMB po- 
larization data. To do this, we adopt and ex t end th e 
TV-point algorithms developed bv lEriksen et al.1 (|2004bf ). 

2. METHODS AND DEFINITIONS 

2.1. N-point correlation functions 

An TV-point correlation function Cm of a stochastic 
field X{n) is defined as 



Cn (8i 



fc ) = (X(n 1 )X(n 2 )---X(n JV ) 



(1) 



where {hi, . . . , njv} spans an TV-polygon defined by k pa- 
rameters, {6*i, . . . , 6k] ■ In the most general case in which 
no assumptions are made concerning the statistical prop- 
erties of the field X, one needs k = 2TV parameters to 



2 



3 N°c 3 



e 



CM 



, a. 



e 



/ou 



Fig. 1. — Center of mass and rotation angles of a triangle. Each 
point i is rotated by an angle on from the global coordinate system 
into the local system defined by the triangle's center of mass. 



uniquely describe such a polygorQ, namely, the individ- 
ual positions of each vertex of the polygon. However, 
in many applications one assumes that X is isotropic 
and homogeneous, and one can therefore average over 
position. In such cases, the number of parameters are 
reduced by three on a two-dimensional surface, corre- 
sponding to translation and rotation. 

In this paper, we will consider only two- and 
three-point correlation functions defined on the two- 
dimensional sphere, and we therefore need, respectively, 
one and three parameters to describe our polygons (ie., 
line and triangle). In the two-point case, the only 
natural parameter choice is the angular distance, 9 — 
arccos(ni -712) between the two points, while in the three- 
point case there is some freedom to choose. For simplic- 
ity, we parameterize the triangle by the lengths of the 
three edges, (#1, #2, #3), where the edges are ordered such 
that B\ corresponds to the longest of the three edges, and 
the remaining are listed according to clockwise traversal 
of the triangle. 

2.2. Polarized correlation functions 

Our goal in this paper is to measure the TV-point cor- 
relation functions of the polarized WMAP maps. To 
do so, we have to ge neralize the methods described by 
lEriksen et al.l (|2004bl) . in order to account for the fact 
that the CMB polarization field is a spin-2 field. Explic- 
itly, the full CMB fluctuation field may be described in 
terms of the three Stokes parameters, l(n), Q(n), and 
U(n). Here I is the usual (scalar) temperature fluctua- 
tion field, and Q and U are two parameters describing 
the linear polarization properties of the radiation in di- 
rection n. 

According to standard CMB conventions, Q and U are 
defined with respect to the local meridian of the spheri- 
cal coordinate system of choice. However, Q and U form 
a spin-2 field, and this means that if one performs a rota- 
tional coordinate system transformation, the Stokes pa- 

12 In this paper, we restrict our interest to fields defined on the 
two-dimensional sphere. 



rameters after transformation becomes 

no \ /i\ 

= cos 2a sin 2a Q , (2) 
\0 -sin 2a cos 2a J \U J 

where a is the local rotation angle between the two 
coordinate systems. iZaldarriaga fc Seliakl ()1997t ) and 
IKamionkowski et al.l ()1997[ ) give complete descriptions of 
the statistics of these fields. 

This property has some important consequences when 
computing TV-point correlation functions: since we as- 
sume that the CMB field is isotropic and homogeneous, 
the value we obtain for a given element of the correla- 
tion function should not depend on the coordinate sys- 
tem. One cannot therefore blindly adopt Q and U as the 
quantities to correlate, since these are in fact coordinate 
dependent. 

For the two-point function, this problem is convention- 
ally solved by defining a local coordinate system in which 
the local meridian passes throu gh the two points of inter- 
est (|Kamio nkowski "et al.l 119971 ). The Stokes parameters 
in this new "radial" system are denoted by Q r and U r , 
and measure the polarization either orthogonally or per- 
pendicularly (Q r ) to the connecting line, or 45° rotated 
with respect to it (U r ). The two-point function defined 
in terms of Q r and U r then becomes coordinate system 
independent. 

We now generalize this idea to higher-order TV-point 
correlation functions. As with the definition of the TV- 
point polygon parameters, we also here have some free- 
dom when choosing the reference points for the coordi- 
nate independent quantities. For instance, two natural 
choices for Q[ and U[ are to define them either with re- 
spect to the edges or with respect to the center of mass 
of the polygon, 

Ei n i fn\ 
"CM = ^ ; | - (3) 

In this paper we choose the latter definition. The result- 
ing geometry is illustrated in Figure [TJ 

Given these new rotationally invariant quantities, X 6 
{I, Q r ,U r }, defined with respect to the local center of 
mass, the polarized correlation functions are simply de- 
fined as 



k ) = (X 1 (n 1 )X 2 (n 2 )---X JV (n J v)) 1 (4) 



Note that in the case of the two-point function, there are 
six independent polarized correlation functions (II, IQ r , 
HT, Q r Q r , Q r U r and U r U r ), while in the three-point 
function case, there are 27 independent components (III, 
IKf, . . ., lTU r Q r , U r U r lT)- 

Algorithmically, we compute these functions from a 
given pixelized sky map, m(h), simply by averaging the 
corresponding products over all available pixel multiplets 
that satisfy the geometrical constraints of the polygon 
under consideration, 



PN(i)' 



. ( 5 ) 

Here i = l,...,TV p runs over the number of available 
pixel multiplets, and Pj(i) is the j'th pixel in the i'th 
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Fig. 2. — Comparison between the analytically calculated two-point function (black) and the two-point function computed from the map 
with the corresponding C;'s (red). The shown functions are, from left to right, the IQ, QQ, and UU correlation functions. 



pixel multiplet. For full details on how to find the rel- 
evant pixel m ultiplets correspondin g to a given polygon 
efficiently, see lEriksen et al.l (|2004b[ ) . 

The correlation functions are binned with a bin size 
tuned to the pixel size of the map. Explicitly, since we use 
a HEALPix resolution of N s n e = 16, with ~ 3° pixels, 
we adopt a bin size of 6°, for a total of 30 bins between 
0° and 180°. 

If all factors in an AT-point multiplet are taken from 
the same map, the resulting function is named an auto- 
correlation; otherwise, it is called a cross-correlation. 
The main advantage of cross-correlations is that the noise 
is typically uncorrelated between maps, while the signal 
(ideally) is strongly correlated. Therefore, the noise con- 
tribution to a cross-correlation averages to zero. For this 
reason, one typically uses auto-correlations as a probe 
of systematic errors (e.g., to check whether the assumed 
noise model is correct), but cross-correlations for a final 
cosmological analysis. In this paper, we consider both 
types. _ 

In Figure [H we show a comparison of three two-point 
correlation functions computed from the same signal- 
only simulation, using two different methods: the red 
curves show the correlation functions computed from the 
power spectrum of the realiza tion, employ ing the ana- 
lytic expression given by, e.g., ISmithl (|2006D ■ The black 
curves show the correlation functions computed directly 
from the pixelized map using the algorithms described 
above. Slight differences are expected due to different 
treatment of pixel windows, but clearly the agreement 
between the two are excellent, giving us confidence that 
our machinery works as expected. 

2.3. Comparison with simulations 

In this paper we perform a standard frequentist anal- 
ysis, in the sense that we compare the results obtained 
from the real data with an ensemble of simulations based 
on some model. Specifically, our null-hypothesis here is 
that the universe is isotropic and homogeneous, and filled 
with Gaussian fluctuations drawn from a ACDM power 
spectrum. Our ensemble consists of N S i m = 100 000 sim- 
ulations. 

For each realization in the ensemble, we compute the 
iV-point correlation functions in precisely the same man- 
ner as for the real data. Finally, we compare the corre- 
lation function derived from the data with the simulated 



functions through a standard % 2 statistic, 

(6) 

b=l ^ 

where the sum runs over all Appoint configurations/bins, 
and the mean fi(b) = (C(b)) and variance cr 2 (b) — 
((C(b) — /i(6)) 2 ) are derived from the simulations. 

However, before computing the y 2 as described above, 
we Gaussianize the distribution of each bin as follows 
(jEriksen et al.ll2005h : 

N sim + 1 J-oc 

where R is the rank of the value under consideration (i.e., 
the number of simulations with a lower value than the 
currently considered), and s is the corresponding Gaus- 
sianized rank. The reason for performing this transfor- 
mation is that the probability distribution for a given 
bin is highly non-Gaussian, and for even-ordered corre- 
lation functions strongly skewed toward high values. A 
direct x 2 evaluation will therefore tend to give too much 
weight to fluctuations that have high values compared 
to fluctuations that have low values. By first explicitly 
Gaussianizing, a symmetric response is guaranteed. 

The final results are quoted as the fraction of simula- 
tions with a higher \ 2 than the real data. By splitting the 
simulations into two disjoint sets, and repeating the % 2 
analysis for each set, we estimate that the Monte Carlo 
uncertainty in the resulting significances is less than 1%. 
For this reason, we quote cases inconsistent with simula- 
tions at more than 99% by instead providing the actual 
number of simulations with higher \ 2 - This allows us 
to distinguish between a case with ~1% significance and 
one with ~ 0.1% significance, but still recognize the im- 
portance of the Monte Carlo uncertainties. Further, we 
conservatively never claim to obtain results with higher 
significance than 99% in any case, despite the fact that 
many results very likely are far more anomalous, as would 
be clear by using more simulations. 

3. DATA AND SIMULATIONS 

In the following, we analyze the five-year WMAP 
sky maps, including both temperature and polarization. 
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TABLE 1 

Statistical significances (two-point functions) 













Correlation 








Frequency 


11 


IQ 


HI 


QI 


QQ 


QU 


01 


UQ 


uu 




KQ85 


KaxKa 










0.40 (0.32) 


0.45 (0.38) 






0.32 (0.25) 


QxQ 


0.08 


0.85 


0.65 




0.01 (348*) 


0.01 (0.01) 






254* (128*) 


QxV 


0.07 


0.26 


0.61 


0.85 


25* (18*) 


0.06 (0.04) 


0.69 


0.01 (338*) 


0.04 (0.02) 


QxKa 




0.81 


0.97 




0.01 (307*) 


0.02 (0.01) 




0.03 (0.02) 


154* (48*) 


QxW 


0.08 


0.84 


0.05 


0.86 


3* (5*) 


0* (0*) 


0.69 


0* (0*) 


253* (222*) 


VxV 


0.07 


0.27 


0.60 




73* (49*) 


0.02 (0.01) 






0.01 (393*) 


VxKa 




0.82 


0.97 




0.01 (357*) 


21* (15*) 




18* (10*) 
37* (8*) 


0.02 (0.02) 
106* (88*) 


VxW 


0.07 


0.84 


0.05 


0.27 


1* (2*) 


42* (14*) 


0.57 


WxKa 




0.78 


0.98 














WxW 


0.07 


0.83 


0.05 




3* (0*) 


71* (34*) 






2* (0*) 
5* (6*) 


WlxWl 


0.09 


0.70 


0.11 




116* (123*) 


0.01 (0.01) 






W2xW2 


0.07 


0.32 


0.22 




95* (83*) 


0.03 (0.03) 






478* (476*) 


W3xW3 


0.07 


0.91 


0.82 




0.78 (0.78) 


0.94 (0.94) 






0.60 (0.60) 


WAxWi 


0.10 


0.53 


0.03 




0* (0*) 


0.08 (0.08) 






1* (3*) 


KQ75 


QxQ 


0.07 


0.96 


0.61 














QxV 


0.06 


0.30 


0.67 


0.97 






0.62 






QxKa 




0.87 


0.95 














QxW 


0.06 


0.91 


0.15 


0.97 






0.64 






VxV 


0.06 


0.31 


0.65 














VxKa 




0.88 


0.96 














VxW 


0.06 


0.92 


0.16 


0.32 






0.62 






WxKa 




0.85 


0.96 














WxW 


0.06 


0.92 


0.15 














WlxWl 


0.08 


0.68 


0.16 














W2xW2 


0.06 


0.33 


0.35 














W3xW3 


0.06 


0.93 


0.94 














W4xW4 


0.10 


0.57 


0.10 















Note. — Statistical significances for the two-point functions as computed from the WMAP data. For the W differencing assemblies, 
the first-year maps were used. For the band data, the five-year co-added maps were used. The Monte Carlo uncertainty is 1%. The 
correlations using both the KQ85 and the KQ75 masks are shown. As the pure polarization correlations are insensitive to which I mask is 
being used, these are only displayed under the KQ85 header. An entry with '. . . ' signifies that some other entry in the table contains the 
same information as that entry, while an entry with ' — ' indicates that the information in that entry was not computed. The listed ratios 
indicate the fraction of simulations with a higher x 2 than the observed data. The entries marked with an '*' signifies that the fraction is 
lower than the Monte Carlo uncertainty, but still larger than zero; consequently, the number of simulations (out of 100,000) with a higher 
X 2 than the observed data are shown instead for these entries. The entries in parenthesis are the significances of the noise-only hypothesis. 



These data are available from LAMBDAF^I. including all 
ancillary data, such as beam transfer functions, noise co- 
variance matrices and foreground templates. The total 
data set spans five frequencies (23, 33, 41, 61, and 94 
GHz), corresponding to the K, Ka, Q, V, and W bands. 

Our main objects of interest are the frequency co- 
added and foreground-reduced sky maps, which are pro- 
vided in the form of HEALPbQ pixelized sky maps. The 
temperature sky maps are given at a HEALPix resolu- 
tion of N s idc = 512, while the polarization maps are given 
at iVgido = 16. To bring the data into a common format, 
we therefore degrade the temperature component to the 
same resolution as the polarization maps, simply by av- 
eraging over sub-pixels in a low-resolution pixel. 

For temperature we consider the Q, V, and W bands, 
and for polarization also the Ka band. This difference 
mirrors the official cosmological WMAP analysis, which 
also uses the Ka band for polarization, but considers the 
same band to be too foreground contaminated in tem- 
perature. 

Our simulations are generated as follows: we first 
draw a random Gaussian CMB realization from the 



13 http://lambda.gsfc.nasa.gov 

14 http://healpix.jpl.nasa.gov 



best-fit five-year WMAP ACDM power spectrum, in- 
cluding multipole moments up to £ max = 1024. This 
realization is then convolved with the instrumental 
beam of each WMAP differencing assembly and the 
iVpi x = 512 HEALPix pixel window, and projected 
onto a HEALPix grid. For temperature, we then add 
a Gaussian noise realization with pixel-dependent vari- 
ance, az — <7q /N i, s (p) , where cXq is the noise variance 
per observation and A^bs is the number of observations 
in that given pixel. The temperature component is then 
degraded to A^dc = 16. 

For the polarization component, we first degrade the 
high-resolution map to A s id = 16, and then add a noise 
term. This is because the polarized noise description 
for WMAP is given by a full noise covariance matrix, 
C p y at A^ido = 16, and not as a simple A^bs count at 
high resolution. The correlated noise realization is gen- 
erated by drawing a vector of standard normal variates, 
r\ ~ AT(0, 1), and multiplying this vector by the Cholesky 
factor of the covariance matrix, n = Ltj, where C = LL*. 

For temperature, we adopt two different masks, namely 
the five-year WMAP KQ85 and KQ75 masks. These ex- 
clude 18% and 28% of the sky, respectively. For the 
polarization component, there is only one relevant mask 
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Angular separation (degrees) Angular separation (degrees) Angular separation (degrees) 

Fig. 3. — Two-point functions computed with the KQ85 mask, as functions of angular separation on the sky. Only the Ka, Q and 
V auto-correlations, as well as the cross-correlations between the W intensity map and the Ka polarization map, are shown. The red 
line shows the two-point function computed from the WMAP five-year co-added maps. The light, medium and dark gray shaded areas 
correspond to la, 2a, and 3a, respectively. The dashed blue line is the maximal point of the correlation histograms, while the black line is 
the theoretically predicted (no noise) two-point function computed from the best-fit five-year WMAP ACDM power spectrum. Each row 
corresponds to a specific combination of modes, while the frequency bands involved are marked on each plot. 
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in the five-year WMAP data release, which excludes 27% 
of the sky. 

Two different polarized foreground templates are used 
in the following. The first is simply the difference be- 
tween the K and Ka bands, smoothed to an effective res- 
olution of 10° to reduce noise, which traces synchrotron 
radiation. The second template is the same starlight tem- 
plate as used for foreground correction in the five-year 
WMAP data release, which traces thermal dust. 

4. RESULTS 

In this section we present the two- and three-point cor- 
relation functions computed from the polarized five-year 
WMAP data. Due to the large number of available and 
unique correlation functions one can form within this 
data set, we plot only a few selected cases, and instead 
show the full set of results only in the form of tabulated 
significances. 

4.1. Two-point correlations 

In Figure [3j we show a selection of two-point corre- 
lation functions derived from the WMAP data. Specif- 
ically, all available Ka, Q and V auto-correlations are 
shown, in addition to the cross-correlations between the 
Ka and W bands. Table Q] shows the \ 2 significances for 
all possible two-point functions from the available data. 
(Note that only unique combinations are actually shown; 
empty table cells indicate that the same combinations 
may be found elsewhere in the table. For example, IQ 
and QI are identical for auto-correlations, but different 
for cross-correlations.) 

Starting with the pure intensity correlations shown in 
the top row of Figure [3l we recognize the by now well- 
known behavior of the CMB temper ature two-point func- 
tion, observed bot h by COBE (e.g.. | Bennett et al.|[l99l) 
and WMAP (e.g., ISpergel"etaLi r2003'l: at small angular 
separations, it is slightly low compared to the ACDM 
model, and at separations larger than 60°, very close to 
zero. This peculiar behavior has attracted the interest 
of both experimentalists and theorists, and some have 
suggested that this could be the signature of a closed 
topological space. However, from the view of a simple 
X 2 test, this function remains consistent with the sim- 
plest ACDM model at the 5%-10% level. 

Next, looking at the temperature-polarization (I-Q/U) 
cross-correlations we see that the overall behavior of the 
Q-, V, and Ka x W band functions are quite different. 
This indicates that the WMAP data still are too noise 
dominated to extract high-sensitivity cosmological infor- 
mation from individual bands. This is also reflected in 
the second and third columns of Table [1] where there is 
a significant scatter between the different results. How- 
ever, we do see that all results are in good agreement 
with the expectations, indicating that the noise model is 
satisfactory. 

The three bottom rows of Figure [3] show the pure po- 
larization correlation functions, which are the main tar- 
get in this paper. And here we see several interesting 
features. First, we recognize a sharp feature in the UU 
correlation function at ~ 141°. This is the angular sep- 
aration between the A- and B sides of the differential 
WMAP detectors, and it was also seen in the pure noise 
temperature two-p oint correlation fun ction in the first- 
year WMAP data (|Eriksen et al.ll2005t) . 



However, even more interesting is the overall very pro- 
nounced large-scale excesses seen in the QQ and UU 
functions: both the Q- and U-band functions lie mostly 
within the 2<7-3<7 confidence regions, and the overall 
agreement with the simulations appears quite poor. 

Again, this is strongly reflected in Table [TJ all QQ 
correlations are anomalous at more than 99% confidence, 
and all UU correlations at more than 98% confidence. In 
the case of the W— band, which happens to be the most 
anomalous of any band, we have good reasons to expect 
such behavior. This band has a sig nificantly higher 1/ f 
knee frequency than any other band (jJarosik et al.f 2003) , 
with the W4 differencing assembly having the highest. 
As a result, the WMAP team has chosen not to use this 
band for cosmological analysis in polarization. But the 
anomalous behavior of the Q and V bands is a priori not 
expected; these are used for cosmological analysis by the 
WMAP team and should be clean. 

We have also generated an ensemble of noise-only sim- 
ulations, and computed the two-point functions from 
these. The corresponding \ 2 fractions are listed in paren- 
theses in Table Q] for the pure polarization modes. Here 
we see that, generally speaking, the noise-only hypothesis 
performs almost, but not quite, as well as the signal-plus- 
noise hypothesis. This is simply due to the fact that the 
WMAP polarization data are strongly dominated by the 
correlated noise, and it is very difficult from a correlation 
function point of view to distinguish between a small sig- 
nal component and a large-scale noise fluctuation. 

4.1.1. Foreground cross-correlations 

The anomalous behavior seen in the Q- and U-band 
two-point functions clearly needs an explanation. And 
typically, when such unexpected behavior is observed, 
one of the the first issues to consider is residual fore- 
grounds. 

To check whether this may be a relevant issue, we first 
simply compute the cross-correlation functions between 
each of the three frequency bands (Ka, Q, and V) and 
the two available foreground templates (synchrotron and 
dust). This is done both for the real WMAP data and 
the simulated ensembles, and the agreement between the 
(pure CMB + noise) simulations and the WMAP data is 
again quoted in terms of a \ 2 fraction. 

The results from these calculations are as follows: for 
the Ka band, we find that the \ 2 significances are 0.51 
and 0.24 for synchrotron and dust, respectively, indicat- 
ing no significant foreground detection in this case. How- 
ever, for the Q band the corresponding numbers are 0.05 
and ~ 0.005, corresponding to correlations statistically 
significant at 2a and ~ 3rx, respectively. For the V band, 
the numbers are 0.05 and 0.02, significant at 2er or more. 

To study these foreground correlations further, we 
compute the two-point functions from the foreground 
templates directly, and fit these to the observed corre- 
lation functions with a single free amplitude for each 
template, A s and ^4d, by minimizing 

Xl = E ( C °bs(6) ~ c **( b ) - A s C s (b) - A d C d (b)) 2^(6, 6') 

66' 

x (C ohs (b') - C sim (b') - A s C s (b') - A d C d (b')) . 

Here, C s ; m (6) is the correlation function mean and 
S sim (6,6') = ((C(6) - Kb))(C(b') ~ fx(b'))) is the co- 
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Amplitude confidence intervals 



Fig. 4. — Template two-point functions fitted to the WMAP data two-point functions. The red line shows the WMAP data, with the 
gray area being the 1 a confidence intervals. The blue line shows the the mean correlation function of the simulations, while the green line 
shows the mean value plus the foreground contribution. 

pared to the observed functions. First, the red curve 
shows the functions derived from the actual WMAP data, 
with gray bands indicating the la uncertainties derived 
from simulations. The blue curve shows the mean cor- 
relation function of the simulations, C s i m (6), and finally, 
the green curve shows the same, but with the foreground 
contribution added in, Cf g (b) = C s i m (6) + A s C s (b) + 
AdCd(b). Clearly, the latter matches the real data far 
better than the pure CMB-plus-noise hypothesis, with 
the best match seen in the UU correlation function. 
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0.0038±0.0013 
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0.0052±0.0010 
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0.0256±0.0042 
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0.0197±0.0047 
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Note. — Mean values and confidence intervals for 
the best-fit amplitudes of the K — Ka and dust template 
two-point functions relative to the WMAP data two- 
point functions. 

variance matrix, both quantities obtained from the noise 
simulations. The indices 6, b' run over all possible pure 
polarization two-point bins (i.e., both angular bins and 
QQ, QU, and UU correlations). 

The resulting best-fit amplitudes from this calculation 
for Ka-, Q- and V bands are shown in Table [2] Here, we 
again see that the .Ka-band amplitudes are generally sig- 
nificantly lower than those of the Q and V bands. (Note 
that the uncertainties quoted in this table only include 
statistical errors, not systematic errors. The significances 
should therefore not be considered as true detection lev- 
els, but are only suggestive.) 

In Figure |4l the fitted correlation functions are com- 



4.2. Three-point correlations 

We now turn to the three-point correlation functions 
derived from the five-year polarized WMAP data. How- 
ever, we note that the three-point function is primarily 
used in the literature as a test of non-Gaussianity. In 
the present case, this will not be case, since we have 
already seen that the data appear to be foreground con- 
taminated, and even the two-point correlation function 
fails a standard % 2 test. The following three-point anal- 
ysis will therefore essentially only be a consistency check 
of the above results, and a demonstration of the general 
procedure, to be further developed in the future, when 
cleaner data become available. 

We consider only cross-correlations of pure polariza- 



TABLE 3 
Three-point function results 



Frequency 



QQQ 



QQU 



QUQ 



Correlation 
QUU UQQ 



UQU 



UUQ 



UUU 



Equilateral functions 



KaxQxV 
KaxVxQ 



0.04 (0.03) 
0.45 (0.40) 



0.31 (0.26) 
0.10 (0.08) 



18* (4*) 
0.07 (0.05) 



0.04 (0.03) 
0.01 (446*) 



0.16 (0.13) 
359* (241*) 



0.01 (416*) 
0.35 (0.31) 



0.07 (0.05) 
0.01 (0.01) 



0.21 (0.17) 
0.11 (0.09) 



Collapsed functions 



KaxQxV 
KaxVxQ 
QxKaxV 
QxVxKa 
VxKaxQ 
VxQxKa 



0.49 (0.45) 
323* (240*) 
0.01 (0.01) 
0.05 (0.04) 
0.21 (0.17) 
0.01 (0.01) 



0.01 (0.01) 
0.01 (0.01) 
0.28 (0.25) 
0.09 (0.07) 
0.04 (0.03) 
0.03 (0.02) 



0.23 (0.18) 

2* (1*) 
0.01 (329*) 
0.01 (0.01) 
180* (5*) 
0.58 (0.51) 



0.03 (0.02) 
0.05 (0.03) 
0.17 (0.14) 
0.10 (0.08) 

0* (2*) 
439* (301*) 



0.63 (0.56) 
0* (1*) 

0.01 (298*) 

0.04 (0.03) 
3* (2*) 

0.01 (0.01) 



0.01 (384*) 
0.11 (0.09) 
231* (179*) 
0.07 (0.06) 
127* (79*) 
0.38 (0.33) 



0.21 (0.18) 
0.05 (0.03) 
0.07 (0.05) 
0.12 (0.09) 

80 (45*) 
0.18 (0.13) 



0.43 (0.38) 
17* (11*) 
0.15 (0.12) 
0.12 (0.09) 
0.01 (456*) 
0.05 (0.04) 



Note. — Statistical significances of the polarization-only three-point functions computed from the WMAP five-year co-added Ka, 
and V maps. The Monte Carlo uncertainty is 1%. See Table [T] for feature explanation. 



tion modes for the cosmologically interesting frequency 
channels. This leaves us with only eight independent cor- 
relation function modes (QQQ, QUQ, UQQ, UUU etc.), 
and only one frequency combination (KaxQxV). Fur- 
ther, we consider only two general types of three-point 
configurations, namely, equilateral and pseudo-collapsed 
triangles. The latter is defined such that two legs of the 
triangle are of equal length, while the third leg spans an 
angular distance of less than 12°. 

The results from these calculations are tabulated in Ta- 
ble [31 and a few arbitrarily chosen modes of the equilat- 
eral three-point functions are plotted in Figure [5l (Plot- 
ting all functions would take too much space, without 
adding any particular new insights.) 

Once again, we see that the x 2 ' s obtained from the 
real data are high, but generally not quite as striking as 
in the two-point case. This is at least partly due to the 
fact that the Ka band is involved in all configurations, 
as well as the fact that the three-point function is always 
"noisier" than the two-point function. We also see that 
the noise-only hypothesis always has higher x 2 's than 
the signal-plus-noise hypothesis, again following the two- 
point behavior. 

It is difficult to interpret these results much further, 
given that there are clear problems already at the two- 
point level. We therefore leave a full interpretation of 
the three-point function to a future publication, when 
cleaner data, either from WMAP or Planck, have been 
made available. 

5. CONCLUSIONS 

We compute for the first time the polarized two- 
and three-point correlation functions from the five-year 
WMAP data, with the initial motivation of constraining 
possible non-Gaussian signals on large scales in the CMB 
polarization sector. However, our main result is a detec- 



tion of a signature consistent with residual foregrounds 
in both the Q- and U-band data, significant at 2a-3a. 

This should not come as a complete surprise, though, 
as a similar result was found in the three-year 
WMAP data t hroug h a direct template fit approach by 
lEriksen et al.1 (j2007f ). They used a Gibbs sampling im- 
plementation to estimate the joint CMB power spectrum 
and foreground template posterior, and found non-zero 
template amplitudes at similar significance levels when 
marginalizing over the CMB spectrum. Given the re- 
sults found in the present paper, this still appears to be 
an important issue for the five-year WMAP data release. 
Resolving this question is of major importance for the 
CMB community, since many cosmological parameters 
depend critically on the large-scale WMAP polarization 
data, most notably the optical depth of reionization, r, 
and the spectral index of scalar perturbations, n s . 

The question of primordial non-Gaussianity in CMB 
polarization data can not be meaningfully addressed as 
long as this issue is unresolved. The three-point correla- 
tion function results shown in this paper therefore mostly 
serve as a demonstration of the approach, and as a cross- 
check on the two-point results. The three-point function 
will obviously become a more important and independent 
quantity in the future, when higher-fidelity polarization 
data become available. 
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APPENDIX 

PRECOMPUTING CONFIGURATION TABLES 



In or der to compute the two-and-three point functions, we have utilized the method introduced by lEriksen et al.l 
(|2004bD . 



http: / / www.eso.org/ science /healpix/ 
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Fig. 5. — Equilateral three-point functions, computed with the KQ85 mask, as functions of angular separation on the sky. Only a 
representative selection is shown. The red line is the three-point function computed from the WMAP five-year co-added maps. The light, 
medium and dark gray shaded areas correspond to la, 2a, and 3a, respectively. The dashed blue line is the maximal point of the correlation 
histograms. Each row corresponds to a specific combination of modes, while the frequency bands involved are marked on each plot. 
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The basic idea of this method is that given the resolution and the mask of the map, and the desired number of 
distance bins, the configurations that correspond to an angle is uniquely determined. One may then do the work of 
finding these configurations only once, saving them in large tables. Each table then corresponds to a certain bin k, and 
the first row in each table contains all the unmasked pixels of the map. Each column contains all pixels at a distance 
given by kd6 from the pixel in the first element of the column. By moving through a column in such a table, one then 
traces out a circle of radius kd0. Thus, by saving the configurations in tables, one pays the initialization costs only 
once. 

An additional, powerful feature of this approach becomes apparent when computing correlation functions of higher 
order than 2. If one thinks of these tables as compasses, tracing out circles of a given length for each pixel, it should be 
clear that, just as one can construct, e.g., triangles using compasses, one can construct triangles on the sphere using 
these tables. If the edges of the desired triangle is of lengths 0, a and /3, one begins with the table, selects a pixel 
from the upper row, and selects a second pixel from the column below this pixel. One then has the baseline of the 
triangle. To construct a triangle with a compass when the baseline is given, one would adjust the compass so that 
it spans a length equal to the length of the second edge of the triangle, and draw a circle with this radius with one 
end of the baseline as the center of the circle. Then, one would again adjust the compass so that it spans the desired 
third length of the triangle, and draw a circle around the other end of the baseline. Then, one would check whether 
the two circles drawn intercept at any point(s). These points would then be the third vertex of the triangle. Carrying 
this analogy to our tables, our two already chosen pixels are the ends of the baseline around which to draw our circles. 
We 'adjust our compass' by finding pixel 1 in the first row of table a, we 'draw a circle around it' by looping over all 
pixels in the column below it, and do the same for pixel 2, expect that we now look in table /3 . By checking whether 
there are any common pixels in these 'circles', we find the triangle(s) whose edges have the desired lengths. 

Since any TV-point polygon can be reduced to triangles, higher-order correlation functions can also be computed 
using this method. 
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